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INTRODUCTION 


The  primary  goal  of  the  1  year  no-cost  extension  to  this  grant  was  to  publish  results.  This  goal  has  been  achieved  with  publication  in 
IEEE  Transactions  in  Medical  Imaging  (1). 

This  project  addressed  the  FY10  PRMRP  subject  of  epilepsy.  Increasing  incidence  of  traumatic  brain  injury  (TBI)  among  soldiers  will 
likely  lead  to  elevated  levels  of  disability  due  to  TBI-related  seizures  and  epilepsy.  The  lack  of  a  reliable  biomarker  hinders  efforts  to 
interrupt  the  evolution  of  epilepsy  from  TBI.  A  new,  network  paradigm  for  analysis  of  brain  imaging  data  suggests  a  new  direction  for 
diagnosing  brain  injury.  Existing  analyses  have  neither  been  applied  to  epilepsy  nor  have  been  validated  by  gold  standard  data.  The 
purpose  of  this  research  was  to  initiate  exploration  of  the  concept  that  network  properties  of  imaging  data  within  the  scope  of 
predicting  the  transition  of  TBI  to  epilepsy.  The  specific  objectives  of  this  project  were  to  develop  a  fast  analysis  protocol  and  validate 
the  analysis  with  gold-standard  invasive  electrophysiology  measurements  from  epilepsy  patients.  The  innovative  aspects  of  the 
research  are  application  of  a  partial  differential  equation  framework  for  fast  analysis,  the  application  of  the  network  paradigm  to 
epilepsy  patients  and  validation  with  gold  standard  invasive  measurements.  The  relevance  of  the  project  to  the  FY 10  PRMRP  topic  of 
epilepsy  stems  from  the  potential  use  of  a  validated,  network  paradigm  as  a  biomarker  of  risk  for  the  transition  from  TBI  to  epilepsy 
with  the  intent  of  interrupting  that  transition. 

BODY 


Research  accomplishments  are  summarized  below  along  the  lines  of  The  Statement  of  Work,  which  outlined  the  following  main  tasks: 

•  Develop  a  partial  differential  equation  (PDE)-based  tractography  methodology  to  enable  fast,  whole-brain  measurements  of 
connectivity. 

•  Validate  noninvasive  measurements  of  connectivity  by  comparison  to  gold  standard,  invasive  electrophysiology 
measurements. 

•  Summarize  and  publish  results. 


Develop  a  partial  differential  equation  (PDE) -based  tractography  methodology  to  enable  fast,  whole-brain  measurements  of 
connectivity. 


The  theoretical  framework  and  tests  of  PDE-based  tractography  have  been 
published  in  a  high  impact  journal  (1).  The  results  go  beyond  a  previously- 
published  conference  paper  (2)  as  follows: 

•  Clarification  of  the  formalism  through  diagrams. 

•  Demonstration  of  efficacy  of  PDE  formalism  when  using  mult- voxel 
target  regions,  as  opposed  to  single-voxel  regions,  setting  the  stage  for 
fast  whole-brain  tractography. 

•  Demonstration  of  PDE  formalism  in  a  wide  range  of  pathways, 
including  corticospinal,  arcuate  and  short  association  fibers. 

•  Demonstration  of  speed  advantage  as  compared  with  Monte  Carlo  and 
graph  theory  based  approaches. 

Additional  work  clearly  demonstrates  the  comparability  of  Monte  Carlo  and 
PDE  approaches  in  the  pathway  connecting  the  posterior  cingulated  and 
entorhinal  cortex  despite  a  12,000  times  acceleration  afforded  by  the  PDE 
approach.  Furthermore,  a  MATLAB  package  has  been  developed  to  implement 
the  tractography. 


Validate  noninvasive  measurements  of  connectivity  by  comparison  to  gold 
standard,  invasive  electrophysiology  measurements. 

PDE  tractography  measures  of  connectivity  were  compared  with 
electrophysiology  measurements  of  cortico-cortico  evoked  potentials  (CCEPS) 
achieved  with  stereotactic  electrodes  (3).  The  use  of  stereotactic  electrodes  avoids  some  of  the  coregistration  error  associated  with  the 
use  of  subdural  grids.  Results  (figure  2)  show  correlation  between  PDE  tractography  measures  of  connectivity  but  not  with  resting-state 
functional  connectivity  (4,  5). 


Figure  1.  Tractography  of  pathway 
connecting  posterior  cingulated  and 
entorhinal  cortex  using  A)  Monte  Carlo  and 
B)  PDE  tractography. 
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Summarize  and  publish  results. 

Work  described  here  has  been  published  in  IEEE  Transactions  on  Medical  Imaging  (1)  and  presented  at  the  2013  Scientific  Meeting  of 
the  International  Society  for  Magnetic  Resonance  (3).  These  are  included  in  the  appendix. 

KEY  RESEARCH  ACCOMPLISHMENTS 

•  Clarification  of  the  theory  of  PDE-based 
tractography. 

•  Development  of  a  MATLAB  package  for 
practical  implementation  of  the  PDE  approach. 

•  Finding  of  correlation  between  gold-standard 
electrode  recordings  and  anatomical  connectivity 
based  on  tractography  while  no  such  correlation 
was  found  by  resting-state  functional 
connectivity. 


REPORTABLE  OUTCOMES 

•  Publication  in  IEEE  Transactions  on  Medical 
Imaging  (1).  Included  in  the  appendix. 

•  Presentation  at  the  2013  Scientific  Meeting  of  the 
International  Society  for  Magnetic  Resonance  in 
Medicine  (3).  Included  in  the  appendix. 


CONCLUSION 

The  primary  goal  of  this  no-cost  extension,  publishing  the  results,  has  been  largely  achieved  with  publication  of  the  method  in  IEEE 
Transactions  of  Medical  Imaging  (1). 

We  have  developed  a  MATLAB  package  that  can  implement  the  PDE  methodology.  This  package  will  be  tested  against  Monte  Carlo 
tractography  and  electrophysiology  in  a  number  of  pathways.  Results  to  date  suggest  that  PDE-based  tractography  measures  of 
anatomical  connectivity  correlate  with  electrophysiological  measures  while  functional  connectivity  does  not.  As  the  speed  of  the 
package  enables  whole-brain  tractography,  the  product  of  this  research  will  enable  network  analyses  of  anatomical  connectivity. 

As  a  scientific  or  medical  product,  the  work  accomplished  represents  a  step  towards  totally  non-invasive  evaluation  of  the  brain  at  risk 
for  epilepsy.  Such  an  evaluation  would  enable  rapid  evaluation  of  pharmacologic  interventions  and  development  of  new  therapies. 
Unfortunately,  although  victims  of  traumatic  brain  injury  are  at  high  risk  for  developing  epilepsy,  there  is  no  clear-cut  way  to  predict  or 
evaluate  strategies  for  treating  epileptic  seizures.  On  the  near  term,  surgical  intervention  for  phamacoresistant  epilepsy  often  relies  on 
highly  invasive  electrode  monitoring  that  is  an  option  for  only  the  most  highly  motivated  patients.  Progress  toward  noninvasive 
detection  of  targets  for  surgical  resection  would  relieve  the  burden  of  suffering  among  these  patients  while  opening  up  new  treatment 
options. 
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Figure  2.  PDE-based  connectivity  (A)  correlates  with  cortico- 
cortico  evoked  potentials  (CCEPs)  while  resting  state  functional 
connectivity  (B)  does  not. 
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Logical  Foundations  and  Fast  Implementation 
of  Probabilistic  Tractography 


Myron  Zhang,  Ken  E.  Sakaie,  and  Stephen  E.  Jones 


Abstract — Although  tractography  can  noninvasively  map 
axonal  pathways,  current  approaches  are  typically  incomplete  or 
computationally  intensive.  Fast,  complete  maps  may  serve  as  a 
useful  clinical  tool  for  assessing  neurological  disorders  stemming 
from  pathological  anatomical  connections  such  as  epilepsy.  We 
re-frame  tractography  in  terms  of  logic  and  conditional  prob¬ 
abilities.  The  formalism  inherently  includes  global  constraints 
and  can  compute  connections  between  any  two  arbitrary  regions 
of  the  brain.  The  formalism  also  lends  itself  to  a  fast  implementa¬ 
tion  using  standard  partial  differential  equation  solvers,  which 
makes  whole-brain  probabilistic  maps  of  anatomical  connectivity 
feasible.  We  demonstrate  results  of  our  implementation  on  in 
vivo  data  and  show  that  it  outperforms  Monte  Carlo  approaches 
in  both  computation  time  and  identification  of  pathways. 

Index  Terms — Diffusion  MRI,  probabilistic  tractography, 
connectome,  diffusion  tensor  imaging,  connectivity  analysis 


I.  Introduction 

Tractography  is  the  only  method  that  can  noninva¬ 
sively  map  axonal  pathways  in  the  brain  [l]-[4].  Current 
methods  are  typically  incomplete  or  computationally  inten¬ 
sive.  For  example,  standard  deterministic  streamline  tracto¬ 
graphy  typically  fails  to  delineate  known  lateral  projections  of 
the  corticospinal  tract  (CST)  [5],  [6]  or  transcallosal  connec¬ 
tions  between  hand  areas  of  the  motor  cortex  [7].  Probabilistic 
tractography  [8],  [9]  can  overcome  such  limitations  by 
generating  multiple  paths  over  a  probability  distribution  but 
can  be  computationally  intensive.  Furthermore,  most  algo¬ 
rithms  use  only  local  information.  That  is,  each  step  along  a 
pathway  is  only  determined  by  information  at  the  position  of 
that  step,  ignoring  knowledge  about  distal  regions. 

We  introduce  a  generic  formalism  for  tractography  based  on 
probability  theory.  The  formalism  inherently  includes  global 
information  about  the  terminations  of  pathways  and  the 
presence  of  boundaries.  The  formalism  also  lends  itself  to  fast 
implementation  by  standard  numerical  solvers  for  partial 
differential  equations  (PDE).  After  a  theoretical  derivation  of 
the  formalism,  we  demonstrate  advantages  in  terms  of  speed 

Copyright  (c)  2010  IEEE.  Personal  use  of  this  material  is  permitted. 
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and  identifying  pathways  that  are  otherwise  difficult  to 
delineate.  The  advantages  are  particularly  evident  when 
assessing  long  pathways,  such  as  transcallosal  connections 
between  hand  motor  regions,  pathways  not  lying  along  major 
tracts,  such  as  the  arcuate  fasciculus  and  pathways  near 
boundaries,  such  as  U-fibers  connecting  cortical  regions  along 
gyri. 

The  proposed  method  uses  probability  theory  differently 
from  previous  constructions.  Standard  formulations  of 
probabilistic  tractography  generate  an  ensemble  of  tracks  [8], 

[9] .  A  connectivity  map  is  generated  by  showing  the  number 
of  tracks  in  each  voxel.  A  given  track  in  the  ensemble  consists 
of  a  sequence  of  line  segments  that  connect  to  form  a  stream¬ 
line.  The  line  segments  initiate  at  a  seed  point,  the  direction  of 
each  segment  is  generated  in  a  probabilistic  fashion,  and  the 
streamline  terminates  upon  encountering  a  stopping  criterion. 
Furthermore,  the  probability  distributions  at  each  voxel  are 
derived  and  sampled  using  only  local  information.  For 
example,  the  direction  of  a  given  segment  is  not  influenced  by 
whether  or  not  subsequent  segments  intersect  a  target  region 
of  interest. 

The  method  proposed  here  focuses  on  the  expected  number 
of  tracks  passing  through  each  voxel,  rather  than  the  trajectory 
of  segments  along  the  track,  from  a  probabilistic  standpoint. 
Tracks  are  not  explicitly  generated,  but  the  number  of  tracks  is 
required  to  be  continuous.  Furthermore,  the  local  information 
from  imaging  data  is  conditioned  by  requiring  all  the  implicit 
tracks  to  intersect  both  a  seed  and  target  region  and  avoid  the 
tissue  boundary.  As  standard  probabilistic  tractography 
generates  populations  of  tracks  by  random  sampling  methods, 
we  will  refer  to  it  as  Monte  Carlo  (MC)  tractography.  The 
method  introduced  here,  implemented  with  PDE  techniques, 
will  be  referred  to  as  PDE  tractography. 

Although  the  PDE  approach  was  designed  to  mimic  the 
behavior  of  track  counts  generated  from  MC  tractography,  the 
resulting  formalism  bears  some  similarity  to  graph-based 
shortest  path  [10],  [11]  and  flow-based  approaches  [12],  [13]. 
Differences  with  regard  to  interpretation  and  performance 
among  the  approaches  will  be  examined  with  direct  compari¬ 
son  of  results  demonstrated  with  the  graph-theoretic  approach 

[10] . 

II.  Theory 

The  goal  is  to  calculate  the  number  of  tracks  in  a  given 
voxel  subject  to  the  conditions  that  the  tracks  leave  a  seed 
region  and  terminate  in  a  target  region  before  intersecting  a 
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Fig.  1.  Illustration  of  notation.  Tracks  passing  through  a  voxel,  x,  pass 
through  a  neighbor  voxel,  x*.  For  clarity,  only  two  neighbors  are  shown.  Solid 
lines  indicate  tracks  originating  in  a  seed  voxel,  s,  and  terminating  at  a  target 
voxel,  t,  without  intersecting  the  boundary.  Dashed  line  indicates  a  track  that 
hits  the  boundary  before  reaching  the  target.  Although  all  tracks  satisfy  the 
condition  XJC  (track  passes  through  neighbor  x,-  before  entering  x),  only  the 
tracks  indicated  by  solid  lines  satisfy  the  condition  XtT  (track  terminates  at 
target). 


boundary.  To  achieve  this  goal,  we  will  first  develop  a  generic 
formalism  for  considering  the  track  number  from  the  stand¬ 
point  of  conditional  probabilities  and  logic.  We  will  then 
relate  this  formalism  to  diffusion  MRI  data  and  a  practical 
algorithm  for  calculating  the  number  of  tracks. 

The  notation  follows  conventions  for  conditional  probabili¬ 
ties  that  have  been  used  in  other  treatments  of  the  tractography 
problem  [8],  [14].  p(A\B,C)  denotes  the  probability  of 
proposition  A  given  the  condition  that  propositions  B  and  C 
are  true.  The  comma  represents  the  logical  AND  operation. 
The  derivation  uses  the  product  rule  [15]: 

p{A,B\c)  =  p{A\B,C)p{B\c).  (1) 

We  will  use  lowercase  letters  to  indicate  voxels  and  uppercase 
letters  to  indicate  conditions.  For  example,  t  indicates  a  target 
voxel,  x  a  voxel  along  a  track  and  x,-  a  neighbor  of  x.  XfT 
indicates  the  proposition  that  a  track  at  voxel  xz  eventually 
reaches  the  target  voxel,  t ,  and  XJC  indicates  the  proposition 
that  a  track  moves  from  voxel  x,-  to  x  (Fig.  1).  In  addition,  we 
define  cp(x)  as  the  number  of  tracks  at  voxel  x. 

Continuity  arguments  lead  to  the  central  equations  for  the 
derivation.  The  number  of  tracks  in  a  voxel  is  directly  related 
to  the  number  of  tracks  in  neighboring  voxels: 

(P{x)  =  YJp{XiX\XiT,l)(P{xi),  (2) 

i 

which  states  that  tracks  in  a  voxel,  x,  arrive  via  a  neighboring 
voxel,  Xi  and  that  only  tracks  terminating  at  the  target,  t ,  are 
considered.  The  relation  is  through  the  transition  probabili¬ 
ties,  p(XtX  |  XT,  I) ,  of  a  track  moving  from  a  neighbor,  xh  to 
the  voxel  x  subject  to  the  condition,  X{T,  that  the  track 
eventually  reaches  the  target.  The  symbol  /  denotes  all  prior 
information,  including  the  implicit  condition  that  the  tracks 
start  at  a  seed  voxel  and  do  not  intersect  a  boundary.  The 
conditions  introduce  global  constraints  by  enforcing  tracks  to 
intersect  the  target.  In  a  typical  derivation,  only  local  transi¬ 
tion  probabilities,  dependent  only  on  properties  at  the  voxel  x 
and  the  position  of  its  neighbor,  xh  would  be  used  without  any 
inherent  conditions  regarding  the  target  or  boundary. 

To  calculate  the  global  transition  probabilities  in  (2),  we 
first  solve  for  p(XT  \  I) .  These  probabilities  are  also  subject  to 


Fig.  2.  Illustration  of  joint  probabilities  in  (5).  Tracks  passing  through  voxel 
Xi  will  pass  through  a  neighbor  voxel,  x,  or  some  other  voxel,  represented  by 
x  .  Tracks  represented  by  both  lines  satisfy  the  condition  XtT.  However,  only 
tracks  represented  by  the  solid  line  satisfy  the  joint  condition  XjXJCT,  which  is 
therefore  equivalent  to  the  joint  condition  XjX,XT. 


a  continuity  condition: 

p{XT\l)  =  Ydp{XXi\l)p{XlT\l),  (3) 

i 


and  involve  the  local  transition  probability  of  moving  from  x 
to  xh  p(XXi  1 1) .  These  local  transition  probabilities  are 
derived  from  the  diffusion  MRI  data.  The  approach  specified 
here  involves  integration  over  a  fiber  orientation  distribution 
function  (fODF)  [16].  The  fODF  is  integrated,  on  a  voxel-by- 
voxel  basis,  over  sub-regions  of  the  unit  sphere  nearest  a  line 
connecting  the  center  of  a  voxel  with  each  neighbor  voxel. 
Details  are  provided  in  the  Methods  section.  If  a  diffusion 
orientation  distribution  function  (dODF)  [17]-[19]  is  used 
instead,  the  resulting  map  describes  the  anisotropic  random 
walk  of  water  molecules. 

Equations  (2)  and  (3)  are  finite  difference  equations  for 
which  a  number  of  standard  and  fast  numerical  solvers  exist. 
Solution  requires  us  to  specify  the  transition  probabilities 
p(XX\XT,I)  to  solve  for  cp(x)  in  (2)  and  p(XXi  | /)  to 
solve  for  p(XT  |  /)  in  (3). 

We  will  now  use  logical  considerations  to  derive  the  global 
transition  probability,  p(XtX  \  XT,  I) .  From  the  product  rule 
(1),  we  find: 


P(X,X\X,TJ) 


p{x,x,x,t\i) 

p(XiT\l) 


(4) 


Physical  considerations  (Fig.  2)  lead  to  the  relation: 

p{XiX,XiT\l)  =  p{XiX,XT\l).  (5) 


Although  proposition  XjX  is  not  independent  from  proposition 
XjT ,  propositions  X^C  and  XT  are  independent  if  we  assume  the 
transition  process  is  Markov,  so: 


p(x,x,xr^)  =  p{xtx\i)p{xr\i)  (6) 


Inserting  (6)  into  (4),  we  find: 


p{X,X\X,T,l) 


p{XiX\l)p{XT\l) 

p(XlT\l) 


(7) 


Given  the  local  transition  probabilities,  p(XX  |  /) ,  one  can 
solve  (3)  to  derive  p(XT  |  /) .  One  can  then  derive  the  global 
transition  probabilities  with  (7)  and  then  solve  for  the  track 
number,  cp(x),  in  (2). 

Essential  to  solving  the  finite  difference  equations  (2)  and 
(3)  is  the  proper  selection  of  boundary  conditions.  For  (2): 

m=o  (8) 

(Pis)  =  N  (9) 
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where  the  number  of  pathways  at  the  boundary,  b ,  is  set  to 
zero  and  number  of  pathways  at  the  seed  voxel,  5,  is  set  to  a 
number  N.  These  conditions  directly  reflect  physical  con¬ 
straints.  N  represents  the  number  of  tracks  originating  at  the 
seed  that  subsequently  reach  the  target  without  hitting  the 
boundary.  N  acts  as  an  overall  scaling  factor  for  the  number  of 
tracks,  <p(x),  found  at  each  voxel  in  the  imaging  volume.  For 
example,  if  we  set  N=  100  and  cp(x)  =  20,  that  means  20%  of 
the  tracks  between  s  and  t  pass  through  x.  In  all  of  our  results, 
we  set  N  =  1,  so  that  cp{x)  represents  the  fraction  of  tracks 
passing  through  x.  We  refer  to  this  scaled  track  number  as 
“track  density.” 

For  (3),  we  have: 

P(BT\I)  =  0  (10) 

P(TT\I)  =  1.  (11) 

The  first  condition  excludes  tracks  at  the  boundary.  The 
second  is  required  by  consistency — tracks  at  the  target  clearly 
reach  the  target.  We  also  impose  the  constraints: 

P(TT  |  /)  =  0  (12) 

P(S,S\I)  =  0  (13) 

for  all  neighbors  tt  of  the  target,  t,  and  neighbors  st  of  the  seed, 
s,  when  solving  (2).  Equation  (12)  reflects  the  fact  that  tracks 
do  not  emanate  from  the  target,  while  (13)  prevents  tracks 
from  returning  to  the  seed  after  they  have  left  the  seed. 

The  final  mathematical  expressions  for  numerical  computa¬ 
tion  of  (2)  and  (3)  are  two  large,  linear  systems  of  finite- 
difference  equations,  which  can  be  written  as  square  matrices 
whose  size  is  determined  by  the  number  of  computational 
voxels  used  to  represent  the  brain  volume.  Since  the  value  of 
any  given  voxel  is  determined  solely  by  the  values  of  adjacent 
voxels,  the  matrices  are  sparse.  The  equations  are  finite- 
difference  forms  of  PDEs  [20],  suggesting  that,  in  the 
continuum  limit,  they  could  be  derived  from  a  PDE. 

III.  Methods 

We  compare  the  PDE  approach  to  a  Monte  Carlo  approach 
and  a  shortest  path  approach  with  in  vivo  data.  Tracking 
performance  is  evaluated  by  comparing  the  similarity  of  track 
density  maps  and  track  geometries  between  the  approaches 
and  by  comparing  computation  time.  We  also  construct  a 
digital  phantom  to  evaluate  discretization  error  and  demon¬ 
strate  key  conceptual  and  performance  points  in  the  PDE 
method. 

A.  In  Vivo  Data 

In  vivo  data  were  acquired  under  a  protocol  approved  by  the 
local  internal  review  board.  Five  subjects  were  imaged  on  a 
Siemens  TIM  Trio  (Siemens  Medical  Solutions,  Erlangen, 
Germany)  with  a  standard  12-channel  head  coil.  The  HARDI 
acquisition  provided  whole-brain  coverage  with  2.5  mm 
isotropic  voxels  (256  mm  x  256  mm  FOV,  102  x  102  matrix, 
48  slices.  TE  -  77  msec,  TR  -  6500  msec,  BW  -  1442 
Hz/pixel,  partial  Fourier  factor  =  5/8,  61  non-collinear 
diffusion-weighting  gradients  with  robust  ordering  [21]  with  b 
=  1000  sec/mm2  and  7  b  =  0  volumes,  2  averages).  Motion 
correction  was  performed  with  an  iterative  algorithm  [22]  that 


updated  gradient  vectors  [23]. 

For  in  vivo  data,  seed  and  target  points  were  placed  to 
examine  several  important  pathways.  Shown  here  are  results 
from  corticospinal  tract,  transcallosal  pathway  between  hand 
regions  of  motor  cortex,  arcuate  fasciculus,  superior  longitudi¬ 
nal  fasciculus,  and  a  short  association  fiber  (U-fiber).  The 
length  of  and  presence  of  multiple  fiber  crossings  along  the 
transcallosal  pathway  between  hand  regions  of  motor  cortex 
make  this  pathway  hard  to  define,  even  by  probabilistic 
tractography.  The  arcuate  fasciculus  constitutes  the  inferior 
margin  of  the  superior  longitudinal  fasciculus.  Each  of  these 
pathways  can  be  challenging  to  identify  by  tractography  due  to 
their  length  and  proximity  to  the  boundary.  The  superior 
longitudinal  fasciculus  and  short  association  fiber  are  used  to 
compare  utility  of  the  PDE  method  in  long  and  short  associa¬ 
tion  pathways. 

In  all  tracts  other  than  the  corticospinal  tract,  a  single-voxel 
seed  and  target  was  used  to  define  the  pathway.  For  the 
corticospinal  tract,  both  single-voxel  and  multi-voxel  regions 
were  used  to  demonstrate  the  extensibility  of  the  PDE  method 
to  multi-voxel  seeds  and  targets.  The  seed  included  569  voxels 
in  the  motor  strip,  segmented  by  FreeSurfer  [24],  and  the 
target  included  4  voxels  in  the  brainstem. 

All  in  vivo  tractography  (PDE,  MC,  and  shortest  path)  was 
performed  with  a  white  matter  mask,  generated  with  SPM 
[25],  to  define  boundaries.  For  the  PDE  tractography  to  run, 
all  that  is  needed,  in  principle,  is  a  tissue  mask  to  define 
boundaries  and  transition  probabilities  defined  within  the 
mask.  In  practice,  however,  implausible  connections  through, 
for  example,  ventricles  or  across  the  midsagittal  fissure  may 
occur. 

B.  Digital  Phantom 

Sampling  of  the  fiber  orientation  distribution  along  26 
directions,  as  opposed  to  a  continuum  of  directions,  will  lead 
to  discretization  error  when  the  fiber  orientation  lies  oblique  to 
the  imaging  grid.  To  examine  the  effect  of  discretization  error, 
we  constructed  a  simple  digital  phantom  simulating  a  single 
straight  white  matter  fiber  in  the  midst  of  gray  matter.  The 
fiber  is  aligned  in  the  axial  plane  at  angles  with  respect  to  the 
right-left  axis  ranging  from  0°  to  45°.  The  signal  profile  of 
white  matter,  as  a  function  of  diffusion  weighting  gradient,  is 
that  of  an  axially  symmetric  diffusion  tensor  with  fractional 
anisotropy  (FA)  of  0.707  and  mean  diffusivity  (MD)  of  0.7  x 
10"3  mm2/sec  with  principal  eigenvector  aligned  along  the 
fiber.  The  signal  profile  of  gray  matter  is  that  of  isotropic 
tissue  (FA  =  0,  MD  =  0.7  x  10"3  mm2/sec).  Due  to  symmetry, 
properties  for  angles  ranging  from  90°  to  45°  will  be  the  same 
as  angles  from  0°  to  45°. 

The  length  of  the  fiber  was  fixed  at  52  voxel-widths,  and 
the  width  and  thickness  at  1  voxel-width.  The  simulated  data 
has  matrix  size  and  diffusion  gradient  profile  similar  to  that  of 
the  in  vivo  data,  with  dimensions  of  103  x  103  x  48  voxels. 
The  white  matter  fiber  was  centered  in  the  24th  axial  slice. 
The  computational  mask  was  one  voxel  inward  from  all  faces 
of  the  data  matrix  (101  x  101  *46  voxels),  which  included  the 
white  matter  fiber  and  all  surrounding  gray  matter.  As 
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multiple  b  =  0  volumes  are  used  to  counteract  the  impact  of 
noise  for  in  vivo  data  and  no  noise  was  injected  in  the 
simulation,  the  gradient  profile  used  in  the  simulations  used 
only  a  single  b  =  0  image  volume. 

In  addition  to  evaluating  discretization  error,  the  0°  phan¬ 
tom  was  also  used  to  illustrate  p(XT  \  I)  from  (3)  and  the 
importance  of  conditioning  the  probabilities  by  solving  (2) 
with  conditioned  and  unconditioned  transition  probabilities. 
Also,  the  0°  phantom  was  constructed  at  different  resolutions, 
ranging  from  26  x  26  x  12  to  208  x  208  x  96,  to  evaluate  how 
the  speed  performance  of  the  PDE  method  scales  with 
resolution. 

C.  Tractography 

For  the  PDE  method,  local  transition  probabilities  were 
calculated  from  the  fiber  orientation  distribution  calculated  in 
each  voxel  of  in  vivo  data.  The  fiber  orientation  distribution 
function  (fODF)  was  calculated  from  the  HARDI  data  by 
spherical  deconvolution  [16]  with  automatically  optimized 
regularization  [26].  The  resulting  fODF  is  represented  as 
coefficients  of  a  subset  of  spherical  harmonics  that  satisfies 
conditions  of  antipodal  symmetry  and  real,  as  opposed  to 
complex,  values  [27].  Spherical  harmonics  up  to  degree  8 
were  used,  leading  to  45  coefficients.  The  calculation  results 
in  a  model  for  the  orientation  of  white  matter  fibers  in  each 
voxel: 

fODF  (0,</>)  =  YJClmYr(dj)  (14) 

l,m 

where  6  is  the  zenith  angle,  ^  is  the  azimuth  angle,  Y™(0,(/)) 
are  the  modified  spherical  harmonics  and  clm  are  the  coeffi¬ 
cients  specific  to  a  given  voxel.  The  transition  probabilities  are 
the  integral  of  the  fODF  over  a  sub-region  of  the  unit  sphere: 

p{XtX\l)  =  ^fODY{e,4i)ded<l>  (15) 

where  the  integration  ranges  over  the  sub-region  closest  to  the 
vector  connecting  the  centers  of  voxels  x  and  xh  We  perform 
this  integral  numerically  by  calculating  values  of  the  fiber 
orientation  distribution  over  2000  evenly  spaced  points  on  the 
unit  sphere  generated  by  an  Archimedian  spiral 
(http://www.math.niu.edu/-rusin/known-math/97/spherefaq) 
and  summing  values  among  points  nearest  the  unit  vector 
pointing  to  each  of  26  near  neighbor  voxels. 

The  finite  difference  equations  (2)  and  (3)  were  solved  by 
numerical  solution  of  a  large  sparse  matrix  equation  using  the 
bicgstab  function  in  MATLAB  (The  Mathworks,  Natick  MA). 
Successive  over-relaxation  was  also  used,  resulting  in 
equivalent  results  but  requiring  substantially  more  time. 

To  demonstrate  the  importance  of  including  global  informa¬ 
tion  via  conditional  probabilities,  track  counts  were  calculated 
in  a  phantom  using  local  transition  probabilities  only. 
Specifically,  in  (2),  the  conditional  probabilities, 
p(XtX  |  XT  J) ,  are  replaced  by  the  local  transition  probabili¬ 
ties,  p(XiX\I). 

To  compare  the  PDE  method  to  a  MC  method,  we  used  an 
existing  MC  algorithm  based  on  fiber  orientation  distributions 
[7]  adapted  from  [28].  One  billion  tracks  were  initiated  from 
the  seed  point  with  direction  chosen  from  the  fiber  orientation 
distribution  by  rejection  sampling.  The  step  length  was  chosen 


as  0.75  times  the  voxel  length  (1.875  mm),  as  in  the  work  of 
Hagmann  et  al.  [28]  and  in  Lowe  et  al.  [7].  No  optimization 
with  regard  to  step  length  was  performed  in  this  work.  Steps 
with  a  bending  angle  of  more  than  90°  were  reflected  across 
the  origin  to  limit  backtracking.  Streamlines  were  eliminated  if 
they  encountered  the  boundary  or  returned  to  the  seed  before 
reaching  the  target.  Computation  was  performed  using  in- 
house  software  written  in  C  and  implemented  on  a  320-cpu 
cluster. 

We  also  demonstrate  results  generated  from  code  gener¬ 
ously  provided  by  Dr.  Iturria-Medina  for  comparison  with 
graph-based  shortest  path  approaches  [10].  The  code  takes  an 
input  of  connectivity  between  each  voxel  and  its  26  neighbors, 
for  which  we  use  the  same  discretized  fODF  values  as  we  do 
for  the  PDE  method.  The  same  mask  was  used  as  for  all 
approaches. 

For  the  PDE  and  MC  methods,  the  entire  calculation  entails 
1)  calculation,  discretization  and  loading  into  memory  of  the 
fODF,  2)  calculation  of  the  track  density  maps  and  3)  saving 
the  maps  to  disk.  Additionally,  in  the  MC  calculation  track 
density  maps  generated  by  each  CPU  of  the  cluster  are  merged 
in  a  separate  step.  Quoted  PDE  and  MC  calculation  times  are 
only  for  step  2,  calculation  of  track  density  maps.  Step  1 
requires  several  minutes,  but  need  only  be  done  once  for  a 
given  data  set.  Step  3  requires  several  seconds.  Merging 
results  for  the  MC  calculation  from  the  cluster  requires  several 
minutes. 

For  the  shortest  path  method,  the  entire  calculation  entails 
1)  calculation,  discretization  and  loading  into  memory  of  the 
fODF,  2)  generating  the  graph  representation  of  the  brain,  3) 
computation  of  shortest  paths  using  Dijkstra’s  algorithm  and 
4)  saving  the  maps  to  disk.  The  reported  shortest  path 
computation  times  are  only  for  step  3,  computation  of  shortest 
paths  using  Dijkstra’s  algorithm. 

IV.  Results 

First,  we  examine  discretization  error,  intermediate  steps  of 
the  theory  and  the  effect  of  conditioning  transition  probabili¬ 
ties  in  a  single-fiber  digital  phantom.  We  then  present  results 
of  PDE,  MC  and  shortest  path  tractography  in  a  variety  of  in 
vivo  pathways.  Lastly,  we  evaluate  the  speed  performance  of 
PDE  and  shortest  path  tractography  in  increasingly  large 
datasets. 

A.  Phantom  Evaluations 

Figs.  3  and  4  show  the  impact  of  discretization  error  on 
track  density  values.  Track  densities  were  calculated  for  the 
digital  phantom  for  each  integer  angle  between  0°  and  45° 
with  respect  to  the  right-left  axis.  Fig.  3  shows  track  density 
maps  for  angles  of  0°,  22°  and  45°,  from  which  we  can  see 
that  the  results  are  qualitatively  similar.  To  make  a  quantita¬ 
tive  comparison,  a  line  was  drawn  between  the  center  of  the 
seed  and  target  voxels  for  each  angle,  and  values  of  track 
density  in  each  voxel  intersecting  the  line  were  then  assigned  a 
position  coordinate  ranging  from  0  (center  of  seed)  to  1 
(center  of  target).  To  facilitate  comparison  with  the  0° 
configuration,  linear  interpolation  was  used  to  determine 
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Fig.  3.  Visual  comparison  of  the  track  density  maps  of  a  white  matter  fiber 
phantom  rotated  at  0°,  22°,  and  45°  with  respect  to  the  right-left  axis  in  the 
axial  plane.  The  intensity  range  is  set  to  [0,  0.2]. 


Fig.  4.  (a)  Track  density  along  simulated  white  matter  fibers  oriented  along 
the  right-left  axis  (black)  and  in  the  axial  plane  at  a  45°  with  respect  to  the 
right-left  axis  (red).  Fiber  is  1  voxel  in  cross  section.  Position  along  the  track 
is  scaled  so  that  the  seed  is  positioned  at  zero  and  the  target  is  positioned  at  1 . 
(b)  Maximum  relative  error  of  track  density  as  a  function  of  angle.  The  size  of 
error  can  be  large,  but  this  may  be  largely  attributed  to  finite  sampling  and 
linear  interpolation  along  a  very  steep,  nonlinear  curve,  as  seen  in  (a). 

values  at  1000  equally-spaced  points  from  0  to  1.  Fig.  4(a) 
shows  the  profile  of  track  densities  for  0°  and  45°.  Qualita¬ 
tively,  the  difference  is  slight.  We  quantify  the  difference  in 
track  density  cp ,  as  defined  in  (2),  by  calculating  the  maximum 
relative  error  at  angle  6 

max  ^<p(t,6)-<p(t,Q)\l  <p(t,0)},  (16) 

where  t  is  the  position  coordinate  and  6  is  the  angle  with 
respect  to  the  right-left  axis.  The  maximum  relative  error  is 
largest  for  6  =  45°  (Fig.  4(b)).  The  size  of  the  error  can  be 
large  and  occurs  where  the  track  density  exhibits  a  steep  slope 
near  the  seed  and  target.  The  error  may  therefore  result  from 
finite  sampling  and  linear  interpolation  of  track  density  values 
along  the  length  of  the  pathway.  The  cross  section  of  the  fiber 
in  the  simulation  shown  is  1  voxel.  However,  fibers  with 
larger  cross  sections  show  similar  trends  of  comparable 
magnitude,  indicating  that  the  discretization  error  is  inherent 
to  the  algorithm,  not  the  fiber  geometry. 

The  0°  phantom  was  also  used  to  demonstrate  the  results  of 
solving  intermediate  step  (3)  for  p(XT  |  /)  (Fig.  5).  The  value 
at  each  voxel  x  represents  the  probability  that  a  track  from  x 
will  intersect  the  target  but  not  the  boundary.  The  intensity  is 
high  near  the  target  but  falls  off  quickly  with  distance  from  the 
target.  The  rapid  falloff  reflects  the  low  likelihood  of  tracks 
reaching  the  target  before  the  boundary,  a  result  of  the 
boundary  (the  6  faces  of  the  imaging  volume)  being  much 
larger  than  the  target  (a  single  voxel).  The  map  resembles  a 
solution  of  Laplace’s  equation  with  boundary  conditions  of  1 
at  the  target  and  0  on  the  faces  of  the  volume.  However, 
anisotropy  in  the  white  matter  fiber  leads  to  an  oblong  shape 
in  contours,  particularly  near  the  target.  The  value  at  the  seed, 
p(ST  |  /) ,  is  only  3.19><10"4.  Only  about  3  out  of  every  10,000 
tracks  will  reach  the  target  from  the  seed  without  intersecting 


Fig.  5.  The  solution  to  (3),  p(XT  \  I),  in  the  0°  white  matter  fiber  phantom. 
The  intensity  range  is  [0,  1].  The  brightest  white  dot  is  the  target,  the  blue  dot 
represents  the  seed  and  the  boundary  sits  at  the  faces  of  the  imaging  volume, 
including  the  edges  of  the  figure.  Values  fall  rapidly  with  distance  from  the 
target.  Oblong  shapes  of  contours  near  the  target  indicate  effects  of  white 
matter  anisotropy. 


Fig.  6.  (a)  Track  density  map  for  white  matter  fiber  phantom  using 

unconditioned  probabilities  in  the  continuity  equation  (2).  The  intensity  range 
is  [0,  0.2].  The  track  density  appears  isotropic  due  to  the  rapid  falloff  in 
density  in  the  absence  of  conditioning.  The  anisotropy  of  the  fiber  can  be 
barely  discerned  to  the  right  of  the  seed,  (b)  Semi-log  plot  of  track  density 
versus  position  along  the  white  matter  fiber  without  conditioning  on  the  target. 
Position  0  is  at  the  seed  while  1  is  at  the  target. 

the  boundary.  Conditioning  on  the  target  in  (2)  limits  consid¬ 
eration  to  only  the  small  subset  of  tracks  that  reach  the  target 
before  hitting  the  boundary. 

Fig.  6  shows  the  results  of  solving  the  continuity  equation 
(2)  using  the  unconditioned  transition  probabilities  using  the 
same  0°  phantom.  Conceptually,  this  is  equivalent  to  starting 
tracks  from  the  seed  and  then  considering  all  tracks  regardless 
of  whether  they  intersect  the  target  or  boundary.  In  the  track 
density  map  (Fig.  6(a)),  we  can  see  anisotropy  along  the  white 
matter  fiber  in  voxels  close  to  the  seed,  but  as  the  tracks  fan 
out,  the  track  density  decreases  rapidly  with  increasing 
distance  from  the  seed  and  fails  to  rise  back  up  at  the  target.  A 
semi-log  plot  of  track  density  along  the  line  from  seed  to 
target  is  shown  in  Fig.  6(b).  The  conditional  probabilities 
effectively  filter  out  the  tracks  that  strike  the  boundary  or  do 
not  reach  the  target,  leading  to  a  more  focused  track  density 
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Fig.  7.  Maximum  intensity  projection  of  the  corticospinal  tract  generated  by  PDE  (a,  b),  Monte  Carlo  (c,  d),  and  shortest  path  (e,  f)  methods.  Intensity  ranges 
were  set  to  (p  e  [0.05,  0.25]  (a,  b)  and  (p  e  [0.02,  1]  (c,  d),  with  corresponding  color  map  shown  at  the  bottom.  Subsequent  figures  will  only  report  intensity  ranges 
for  q),  with  the  same  color  map  implied. 


TABLE  I 

Relative  Performance  of  PDE,  Shortest  Path  (SP),  and  Monte  Carlo 
(MC)  Tractography, 


Pathway 

CPU- 

seconds 

(PDE) 

CPU- 

seconds 

(SP) 

CPU- 

seconds 

(MC) 

Track 

Efficiency 

(MC) 

Arcuate 

1.4 

8.4 

4.3xio6 

1.4xl0*8 

Transcallosal 

1.1 

8.7 

3.6x10s 

l.lxlO'7 

Longitudinal 

Fasciculus 

0.89 

7.7 

4.3xl06 

1.0x10* 

Short  U-fiber 

0.88 

7.4 

4.4x10s 

6.8x10-“ 

Corticospinal 

(single-voxel) 

1.2 

9.1 

4.6x10s 

2.7x1 0‘6 

Corticospinal 

(multi-voxel) 

1.1 

N/A 

1.6x10s 

2.1X10’5 

For  Monte  Carlo  results,  109  tracks  were  generated  at  the  seed  for  each 
pathway,  and  track  efficiency  denotes  the  fraction  of  those  tracks  that 
terminated  at  the  target.  Comparison  with  multi-voxel  seed  and  target  was 
only  performed  for  PDE  and  MC  tractography  on  the  corticospinal  tract. 

profile  that  can  persist  over  long  distances  (Figs.  3  and  4). 

B.  Anatomical  Results 

Table  I  summarizes  the  computation  times  of  the  PDE,  MC, 
and  shortest  path  methods  in  the  examined  pathways.  We  can 
clearly  see  the  performance  advantage  of  PDE  tractography  as 
compared  with  MC  tractography.  One  billion  tracks  were 
seeded  for  each  pathway  in  MC  tractography.  The  tracking 
efficiency,  or  fraction  of  tracks  leaving  the  seed  voxel  that 
eventually  reach  the  target  voxel,  is  remarkably  low  for  MC 
tractography,  requiring  use  of  a  large  computing  cluster  to 
gain  any  number  of  tracks,  particularly  for  pathways  near  the 
boundary  between  gray  and  white  matter  (arcuate  fasciculus) 
and  for  long  pathways  with  numerous  fiber  crossings  (tran¬ 
scallosal  motor  pathway).  The  computation  time  is  concomi¬ 
tantly  large,  on  the  order  of  1  CPU-month.  While  the  metric  of 


Fig.  8.  PDE  (a,  b)  and  Monte  Carlo  (c,  d)  tracking  between  multi-voxel  seed 
and  target  regions.  Maximum  intensity  projections  are  shown  of  tracking 
between  a  seed  ROI  in  the  motor  strip  (569  voxels)  and  target  ROI  in  the 
brainstem  (4  voxels).  Intensity  range  is  set  to  [0.01,  1]  for  all  images. 

tracking  efficiency  is  not  directly  applicable  to  the  PDE 
approach,  run  time  is  uniformly  on  the  order  of  1  second  for 
any  given  pathway,  in  a  computational  brain  volume  of  65,000 
voxels.  The  shortest  path  approach  requires  less  than  10 
seconds  in  the  same  brain  volume  for  any  pathway. 

Fig.  7  compares  PDE,  MC,  and  shortest  path  tractography 
within  a  readily  defined  pathway,  the  corticospinal  tract 
connecting  brainstem  and  the  hand  knob  of  the  right  motor 
cortex.  Shown  are  maximum  intensity  projections  of  the  track 
density  from  the  PDE  and  MC  methods  and  the  maximum 
likelihood  path  from  the  shortest  path  method.  This  pathway 
typically  poses  few  difficulties  for  tractography,  and  we  find 
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Fig.  9.  Maximum  intensity  projection  of  the  arcuate  fasciculus  generated  by  PDE  (a,  b),  Monte  Carlo  (c,  d),  and  shortest  path  (e,  f)  methods.  Intensity  ranges 
were  set  to  [0.05,  0.3]  (a,  b)  and  [0.01,  1]  (c,  d). 


Fig.  10.  Comparison  of  maximum  intensity  projections  of  the  corticospinal  tract  with  Monte  Carlo  seeding  106  (a),  107  (b),  108  (d),  109  (d)  tracks,  and  the  PDE 
result  (e).  Intensity  ranges  were  set  to  [0.02,  1]  (a-d)  and  [0.05,  0.25]  (e).  As  more  tracks  are  seeded,  the  Monte  Carlo  results  appear  smoother,  approaching  the 
PDE  result. 


Fig.  11.  Maximum  intensity  projection  of  the  transcallosal  pathway  generated  by  PDE  (a,  b),  Monte  Carlo  (c,  d),  and  shortest  path  (e,  f)  methods.  Intensity 
ranges  were  set  to  [0.05,  0.7]  (a,  b)  and  [0.01,  1]  (c,  d).  Monte  Carlo  and  shortest  path  methods  show  a  false  connection  through  the  pons,  which  is  absent  in  the 
PDE  result. 
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that  the  resulting  track  shapes  from  PDE  (Fig.  7(a),  (b)),  MC 
(Fig.  7(c),  (d))  and  shortest  path  (Fig.  7(e),  (f))  methods  are 
similar. 

The  single-voxel  seed  and  target  defining  the  corticospinal 
tract  above  can  be  replaced  by  multiple-voxel  regions  without 
loss  of  generality  or  performance  in  the  PDE  approach.  For 
example,  Fig.  8  shows  the  corticospinal  tract  connecting  a 
seed  region  with  569  voxels  in  motor  cortex  and  4  voxels  in 
brainstem.  The  PDE  computation  time  was  the  same  as  that 
required  for  the  seed  and  target  consisting  of  single  voxels 
(Table  I).  The  track  density  map  shows  details  of  intersection 
with  motor  cortex  (Fig.  8(a),  (c))  that  are  not  apparent  in  the 
single- voxel  seed/target  case  (Fig.  7(a),  (c)).  The  PDE  and  MC 
results  appear  similar,  but  the  PDE  result  shows  more  lateral 
projections  of  the  corticospinal  tract  which  may  explain  the 
larger  extent  of  the  PDE  track  density  map  relative  to  that  of 
the  MC  map. 

Fig.  9  shows  the  arcuate  fasciculus,  which  can  be  difficult 
to  identify  because  it  lies  adjacent  to  gray  matter.  We  can  see 
that  the  PDE  and  shortest  path  results  exhibit  a  similar 
geometry  in  the  axial  view  (Fig.  9(a),  (e)),  but  in  the  sagittal 
view  (Fig.  9(b),  (f)),  the  PDE  track  curves  through  axial  planes 
while  the  shortest  path  track  is  flat,  mostly  restricted  to  a 
single  axial  slice.  The  MC  results  (Fig.  9(c),  (d))  show  a 
similar  geometry  to  the  PDE  results  (Fig.  9(a),  (b)),  but  in 
pathways  close  to  the  boundary  like  this  one,  the  majority  of 
tracks  generated  by  the  MC  method  terminate  prematurely 
upon  intersecting  gray  matter,  as  seen  by  the  patchy  appear¬ 
ance  of  the  track  counts  and  the  low  tracking  efficiency  (Table 
I).  Only  14  tracks  out  of  a  billion  seeded  terminate  at  the 
target. 

Qualitatively,  the  results  of  the  PDE  approach  are  smoother 
than  those  from  the  MC  approach.  This  smoothness  may  arise 
from  two  factors.  First,  the  calculation  of  local  transition 
probabilities  discretizes  the  fODF,  integrating  over  directions 
associated  with  each  of  26  neighbors  of  a  given  voxel.  Such 
integration  smears  the  fODF  with  respect  to  angle,  leading  to 
smoothing  of  resulting  track  density  maps.  Second,  the  PDE 
approach  resembles  the  limit  of  the  MC  approach  as  the 
number  of  samples  approaches  infinity.  It  is  difficult  to 
disentangle  these  two  causes  in  general.  However,  we  can 
demonstrate  how  the  MC  results  approach  the  smooth  PDE 
results  as  the  number  of  MC  samples  increases.  Fig.  10  shows, 
in  the  corticospinal  tract,  how  the  patchiness  of  the  track 
density  gives  way  to  a  smooth  profile  as  the  number  of 
samples  increases.  As  the  MC  approach  does  not  discretize  the 
fODF,  the  smoothing  in  this  case  results  exclusively  from  the 
approach  to  infinity.  We  used  the  corticospinal  tract  rather 
than  the  arcuate  fasciculus  for  Fig.  10  because  the  low 
tracking  efficiency  of  the  arcuate  requires  impractically  long 
computation  times  to  generate  a  sufficient  number  of  tracks. 

Fig.  11  shows  advantages  of  PDE  tractography  in  delineat¬ 
ing  long  pathways  encompassing  multiple  crossings,  such  as 
the  transcallosal  motor  pathway.  The  difficulty  of  MC 
tractography  in  sampling  the  pathway  is  again  shown  by  the 
patchy  appearance  of  the  track  counts  (Fig.  11(c),  (d))  as 
compared  to  PDE  tractography  results  (Fig.  11(a),  (b)). 


Furthermore,  the  MC  tractography  generates  a  well-known  but 
artifactual  “connection”  through  the  pons,  which  does  not 
have  high  track  density  in  the  PDE  result.  This  path  through 
the  pons  is  also  the  maximum-likelihood  connection  generated 
by  the  shortest  path  approach  (Fig.  1 1(e),  (f)). 

To  evaluate  the  PDE  method  in  pathways  of  different 
length,  we  performed  tractography  in  long  and  short  associa¬ 
tion  fibers.  Fig.  A1  (presented  in  the  Appendix)  shows  the 
results  of  tracking  the  superior  longitudinal  fasciculus.  As  in 
the  transcallosal  pathway,  MC  and  PDE  tractography  show 
differences.  As  can  be  seen  in  Fig.  Al,  the  PDE  method 
traverses  the  superior  longitudinal  fasiculus.  For  the  same  seed 
and  target,  the  MC  method  descends  into  the  inferior  longitu¬ 
dinal  fasiculus  and  uncinate  fasciculus,  connecting  with  the 
superior  longitudinal  fasiculus  at  the  ends.  The  shortest  path 
method  track  shows  the  same  geometry  as  the  MC  track. 

Fig.  A2  shows  the  results  of  tracking  a  short  U-fiber.  Due  to 
the  discretization  inherent  to  the  method,  pathways  that  are 
short  compared  to  the  voxel  dimension  may  not  be  as  clearly 
delineated.  This  can  be  seen  from  the  difference  in  appearance 
between  the  long  (Fig.  Al(a),  (b))  and  short  (Fig.  A2(a),  (b)) 
tracks.  While  the  features  of  the  superior  longitudinal 
fasciculus  are  clear,  the  shape  of  the  U-fiber  is  harder  to 
discern  in  the  PDE,  MC  and  shortest  path  methods,  which  all 
present  results  in  discrete  voxels. 

C.  Speed  Evaluation 

Spatial  discretization  inherent  to  the  proposed  method 
imposes  a  strong  limitation  on  spatial  resolution.  Streamline 
tractography  methods  can,  in  principle,  achieve  sub-voxel 
resolution  by  inferring  a  continuous  vector  field  of  fiber 
orientations  from  the  discretely  sampled  field  [4].  Such  a 
result  may  be  approximated  by  the  PDE  approach  by  increas¬ 
ing  voxel  resolution  at  the  acquisition  level  or  by  interpolation 
but  with  a  cost  in  computation  time.  To  assess  the  impact  of 
such  interpolation,  we  constructed  the  0°  phantom  at  eight 
different  resolutions  (26r  x  26r  x  12 r,  r  =  1,  2,  3,  ...,  8).  The 
52  x  52  x  24  voxel  phantom  had  a  mask  containing  6.0x1 04 
tissue  voxels,  approximately  equal  to  that  of  the  in  vivo  tissue 
mask.  The  computation  time  depends  on  the  number  of  tissue 
voxels  within  the  mask,  not  the  total  number  of  voxels  in  the 
imaging  volume.  Fig.  A3  shows  plots  of  PDE  and  graph-based 
computation  times  against  number  of  tissue  voxels  in  the 
mask.  The  run  time  for  the  PDE  method  scales  linearly  with 
the  number  of  voxels  in  the  mask,  while  the  run  time  for  the 
graph-based  tracking  rises  quadratically.  Times  for  the 
phantom  with  a  mask  size  approximately  equal  to  that  of  the  in 
vivo  mask  are  also  indicated.  A  fourfold  increase  in  spatial 
resolution,  corresponding  to  interpolating  in  vivo  voxel 
dimensions  from  2.5  mm  isotropic  to  0.625  mm  isotropic, 
would  increase  the  PDE  computation  time  from  approximately 
1  second  to  4  minutes. 

V.  Discussion 

We  have  introduced  a  theoretical  formalism  for  tractogra¬ 
phy  that  inherently  accounts  for  global  properties  and  lends 
itself  to  fast  implementation.  In  vivo  examples  demonstrate 
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advantages  in  the  ability  to  delineate  pathways  that  are  long 
and  that  course  near  boundaries.  Furthermore,  connections  of 
any  two  points  in  the  brain  can  be  quantified.  While  MC 
tractography  may  yield  zero  streamlines  between  two  distal 
points,  one  may  be  left  to  wonder  if  the  absence  of  streamlines 
is  due  to  an  absence  of  connectivity  or  a  failure  to  generate  a 
large  enough  ensemble  of  streamlines.  The  PDE  approach  will 
always  give  a  nonzero,  but  sometimes  small,  density  of  tracks 
between  two  points.  The  MC  methods  and  PDE  methods 
essentially  solve  the  same  underlying  problem,  but  with 
different  approaches.  Intuitively  the  increased  efficiency  of  the 
PDE  approach  could  be  appreciated  if  the  MC  solution  is 
viewed  as  a  “serial”  computation,  testing  each  possible  track 
one  after  another,  whereas  the  PDE  solution  is  viewed  as  a 
“parallel”  computation,  computing  all  the  possible  pathways 
simultaneously. 

The  use  of  probability  theory  introduced  here  is  inher¬ 
ently  different  from  previously  introduced  approaches  to 
probabilistic  tractography.  Our  approach  is  to  directly 
calculate  the  number  of  pathways  in  a  given  voxel  using 
considerations  of  continuity.  Seminal  work  by  Parker  [9]  and 
Behrens  [8]  use  probabilistic  approaches  to  govern  the 
trajectories  of  each  of  an  ensemble  of  pathways.  The  probabil¬ 
istic  aspect  enters  because  the  step  direction  along  each 
segment  of  each  streamline  is  generated  by  sampling  over  a 
probability  distribution  generated  from  the  diffusion  MRI 
data.  Standard  streamline  tractography  can,  within  this 
framework,  be  seen  as  generating  a  maximum  likelihood 
streamline  from  the  ensemble.  Friman  introduced  a  formalism 
for  calculating  the  posterior  probability  of  individual  stream¬ 
lines  [14].  Similarly,  previously-introduced  global  methods 
[29],  [30]  focus  on  generating  individual  streamlines.  The 
focus  on  direct  calculation  of  track  counts  rather  than  on 
generating  individual  pathways  lends  itself  to  fast  implementa¬ 
tion  by  numerical  PDE  methods. 

The  manner  in  which  the  formalism  incorporates  global 
information  bears  some  discussion.  In  typical,  local  tractogra¬ 
phy,  each  step  along  a  streamline  is  determined  using  only 
information  at  the  location  of  that  step,  ignoring  information 
from  distal  points  along  the  streamline.  In  global  tractography, 
information  from  distal  regions  has  an  influence.  Consider,  for 
example,  a  straight  pathway  in  which  the  fiber  orientation 
halfway  along  the  path  is  corrupted  by  noise,  leading  to  a  fiber 
orientation  perpendicular  to  that  of  other  points  along  the  path. 
In  local  tractography,  limits  on  bending  angle  will  terminate 
tracks  when  they  encounter  the  noisy  voxel.  In  global 
tractography,  the  local  fiber  orientation  information  at  the 
noisy  voxel  is  ignored  in  the  construction  of  a  track  as  long  as 
there  is  consistency  between  the  fiber  orientation  and  the  track 
direction  overall  [30].  In  other  words,  the  global  approach 
departs  from  the  greedy  dominance  of  local  transition 
probabilities  between  neighboring  voxels.  In  the  formalism 
presented  here,  the  transition  probabilities  p(XiX\XiT,I) 
are  constructed  from  local  probabilities,  p(XiX  |  /) ,  that  are 
tempered  by  the  constraints  that  tracks  intersect  the  target  and 
satisfy  the  boundary  conditions,  represented  by  p(XT  |  /)  in 
(7).  The  values  of  p(XT  \  I)  furthermore  encode  what  will 


happen  at  distal  voxels.  For  example,  suppose  that  voxel  xt  has 
large  local  transition  probability,  p(XtX  |  /) ,  but  successive 
transitions  from  x  lead  tracks  away  from  the  target  or  towards 
a  boundary.  Then  p(XT  |  /)  will  be  small,  leading  to  a  small 
value  for  p(XtX  \  XtTJ)  despite  the  large  local  transition 
probability. 

Although  (2)  and  (3)  take  the  form  of  finite-difference 
formulations  of  PDEs,  the  equations  do  not  derive  from  PDEs. 
It  is  of  some  interest,  however,  to  consider  which  PDEs  (2) 
and  (3)  might  approximate.  One  interpretation  of  (3)  is  the 
finite-difference  form  of  Laplace’s  equation  in  a  warped 
space.  To  illustrate,  consider,  in  two  dimensions,  Laplace’s 
equation, 

'  /(x  +  Ay,  y)  -  2  fjx,  y)  +  /(x  -  Ax,  y) 

Jxx  Jyy  ~  2  Ax: 

|  /(x,7  +  Aj)-2/(x,7)  +  /(x,7-A>>) 

2  Ay 

Rearranging  (17)  to  isolate  /(x,y)  gives 

f(x,  y)  » [(A y)f(x  +  Ax,  y)  +  (A y)f(x  -  A x,y)  + 

{Ax) f{x,  y  +  Ay)  +  (A x)f(x,  y  -  A>0]  /  [2{Ax  +  Ay)],  ( 
which  can  be  re-written  in  a  form  analogous  to  (3) 

/(r)  =  Zc(r’r;)/(r;)  (19) 

i 

where  r  corresponds  to  position  (x,y)  and  neighbors  rt 
corresponds  to  neighbors  (x±Ax,y)  and  (x,y±Ay).  The 
coefficients  c(r,r.)  correspond  to  a  function  of  the  grid 
spacing,  either  Ax  /  [2(Ax  +  Ay)]  or  Ay  /  [2(Ax  +  Ay)] ,  that 
varies  from  voxel  to  voxel.  Other  neighbors  can  be  included 
by  considering  higher  order  terms,  and  the  extension  to  three 
dimensions  is  straightforward.  When  Ax  =  Ay  ,  all  coefficients 
are  1/4,  and  the  finite  difference  equation  directly  represents 
Laplace’s  equation.  If  the  coefficients  are  different  for 
different  neighbors  and  at  different  voxels,  the  PDE  is 
equivalent  to  Laplace’s  equation  in  a  warped  metric  space. 
This  interpretation  resembles  ideas  suggested  by  others  [31]- 
[33]  in  different  contexts.  Note  that  the  coefficients  c(r,r.) 
exhibit  antipodal  symmetry  (opposite  directions  are  weighted 
by  the  same  coefficients).  This  symmetry  is  consistent  with 
the  coefficients  p(XXt  \  I)  in  (3),  which  are  derived  from 
antipodally  symmetric  spherical  harmonics.  The  coefficients 
in  (2),  however,  do  not  necessarily  exhibit  antipodal  symme¬ 
try,  so  the  warped-space  interpretation  does  not  apply.  We  do 
not  know  of  any  commonly  used  PDE  that  equivalent  to  (2). 

One  implicit  parameter  of  the  proposed  method  is  the 
number  of  neighbors  defined  for  each  voxel.  On  a  three 
dimensional  grid,  each  voxel  has  6  nearest  neighbors  and  26 
adjacent  neighbors.  The  results  presented  use  26  neighbors. 
Use  of  26  neighbors  may  raise  the  concern  that  the  continuity 
of  track  number  is  violated  because,  when  considering 
transitions  among  26  neighbors,  a  track  may  traverse  more 
than  one  neighbor  voxel,  leading  to  over-counting.  To 
illustrate,  consider  the  simplified  2-dimensional  case  (Fig.  12) 
in  which  a  voxel  x  has  nearest  neighbor  X]  and  an  adjacent,  but 
not  nearest,  neighbor  x2.  If  we  consider  just  the  number  of 
tracks  in  each  voxel,  continuity  is  not  satisfied.  Voxel  x 
contains  three  tracks,  but  the  neighbors  contain  a  total  of  four 
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Fig.  12.  Illustration  of  continuity  when  including  non-nearest  neighbors.  As 
explained  in  the  text,  conditional  probabilities  avoid  over-counting. 

tracks,  two  in  x\  and  two  in  x2.  However,  weighting  with  the 
conditional  probabilities,  p(XiX\XiT,I)  in  the  continuity 
equation  (2)  assures  that  each  track  is  counted  only  in  the 
neighbor  x,-  when  the  track  directly  enters  x  from  xz.  In  the 
example,  all  tracks  in  x\  enter  x  directly,  giving 
p{XxX  |  XxTJ)  =  1 ,  while  only  half  of  the  tracks  in  x2  enter  x 
directly,  giving  p(X2X  \  X2T,I )  =  1  /  2  ,  resulting  in 

<p(X)  =  p(XlX\XlT,l)<p(Xl)  +  p(X2X\X2T,l)<p{X2) 

=  (1)(2)  +  (1  /  2)(2)  =  3, 
satisfying  continuity. 

A.  Adaptability  and  Extensions 

The  methods  introduced  can  be  readily  adapted  to  incorpo¬ 
rate  aspects  of  nearly  any  existing  methodology.  For  example, 
the  orientation  probability  distributions  generated  by  FSL’s 
bedpostx  program  [5],  [34]  can  be  used  to  generate  local 
transition  probabilities  instead  of  fODFs  generated  by 
spherical  deconvolution.  Alternatively,  counts  of  streamlines 
generated  by  diffusion  tensor-based  tractography  can  also  be 
used  [10]. 

The  formalism  may  be  extended  or  modified.  Besides  using 
different  transition  probabilities,  the  boundary  conditions  can 
be  altered  to  specify  the  derivative  of  the  track  count  (Neu¬ 
mann  boundary  conditions)  instead  of  the  value  (Dirichlet)  or 
a  mixture  of  the  two.  Specifying  the  derivative  may  be  more 
appropriate  at  boundaries  where  the  tracks  are  expected  to 
align  parallel  to  the  boundary,  such  as  ventricles  or  brainstem, 
as  opposed  to  boundaries  where  the  tracks  are  expected  to 
terminate  at  the  boundary,  such  as  the  neocortex. 

Our  focus  on  track  counts  is  largely  driven  by  our  long-term 
strategy  of  defining  pathway-dependent  statistics  for  assessing 
diseased  tissue.  For  example,  we  have  previously  assessed 
anatomical  connectivity  along  the  motor  pathway  by  integrat¬ 
ing  metrics  of  tissue  integrity,  such  as  radial  diffusivity,  along 
the  pathway  [7].  For  such  a  strategy,  the  parameterized  path  of 
individual  streamlines  is  not  absolutely  necessary,  but  the 
number  of  streamlines  in  each  voxel  can  be  used  to  generate  a 
weighted  mean  of  tissue  integrity  along  the  pathway. 

However,  it  is  possible  to  construct  streamlines  from  the 
track  density  maps.  The  weighted  sum  of  vectors,  di ,  pointing 
from  the  center  of  a  voxel,  x,  to  the  centers  of  each  of  its 
neighbors,  xh  yields  a  vector: 

A*)  =  Zp(VV|V7VH,  (20) 


a  b 


Fig.  13.  Streamlines  generated  from  the  motor  strip  to  brainstem  using  the 
PDE  formalism.  For  display,  the  3 -dimensional  streamlines  are  projected  onto 
a  coronal  (a)  and  sagittal  (b)  plane. 

where  the  weighting  is  the  conditional  probability  derived  in 
(7).  The  vector  d(x)  can  be  interpreted  as  the  expectation 
value  of  the  direction  of  streamlines.  Defining  d(x)  at  each 
voxel  results  in  a  vector  field  that  can  be  integrated  from 
points  in  the  seed  to  generate  streamlines.  The  resulting 
streamlines  inherit  the  conditions  inherent  to  the  track  density 
map,  such  as  termination  at  the  target  and  avoidance  of 
boundaries.  An  example  of  streamlines  for  the  corticospinal 
tract,  using  the  same  seed  and  target  regions  as  Fig.  8,  is 
shown  in  Fig.  13.  Preliminary  work  suggests  that  these 
streamlines  avoid  unphysical  arrangements  such  as  tight  bends 
and  multiple  visits  to  the  same  voxel.  Future  work  will 
examine  such  streamlines  in  more  detail. 

B.  Comparison  with  Other  PDE  and  Graph  Methods 

A  comparison  in  terms  of  interpretation  can  now  be  made 
with  graph-based  shortest  path  [10],  [11]  and  flow-based 
approaches  [12],  [13].  In  our  PDE  formalism,  the  resulting 
track  count,  tp(x),  can  be  interpreted  as  the  number  of  tracks  in 
a  given  voxel,  x,  that  start  at  a  seed  region  and  terminate  at  a 
target  region  without  intersecting  a  boundary.  The  interpreta¬ 
tion  is  similar  to  connectivity  values  generated  by  the 
probtrackx  program  of  FSL  [8]  if  the  targetmasks"  option  is 
used  and  is  the  default  of  the  MC  algorithm  proposed  by  Lowe 
et  al.  [7]. 

The  interpretation  differs  from  those  of  graph-based  shortest 
path  methods  and  flow-based  algorithms.  Graph-based 
shortest  path  methods  [10],  [11]  construct  a  network  of  edges 
between  each  pair  of  neighboring  voxels.  The  edge  weights 
are  similar,  in  principle,  to  the  local  transition  probabilities 
introduced  here.  The  product  of  edge  weights  along  a  path 
defines  the  weight  of  the  path,  and  the  algorithms  find  the  path 
with  maximum  weight  to  define  connectivity  between  a  seed 
voxel  and  every  other  voxel.  The  result  differs  from  the  PDE 
approach  because  the  edge  weights  are  not  conditioned  on 
termination  at  the  target  and  avoidance  of  boundaries.  As  only 
a  single  path  is  found  between  seed  and  each  voxel,  the  graph- 
based  methods  represents  a  map  of  maximum  likelihood  paths, 
rather  than  an  ensemble  of  pathways. 

Flow-based  approaches  model  fluid-like  movement  through 
a  vector  field  derived  from  the  diffusion  data.  The  resulting 
maps  of  connectivity  are  conceived  as  steady-state  flow 
between  a  source  and  sink  [12]  or  the  arrival  time  of  a  wave 
front  [13].  As  the  approach  of  Hageman  et  al.  finds  a  steady 
state  that  is  consistent  throughout  the  entire  vector  field, 
including  target  and  boundary,  it  bears  strong  similarity  to  the 
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approach  introduced  here.  However,  the  underlying  equations 
differ.  The  approach  of  Hageman  et  al.  [12]  uses  the  Navier- 
Stokes  equations,  literally  modeling  water  movement  through 
the  brain.  The  approach  proposed  here  considers  conservation 
of  track  number  as  the  fundamental  idea.  Front  propagation 
methods  do  not  impose  a  constraint  on  reaching  a  target.  Front 
propagation  methods  can,  however,  incorporate  some  physical 
constraints  in  a  manner  similar  to  streamline  tractography, 
such  as  stopping  propagation  in  areas  of  low  fractional 
anisotropy. 

An  apparently  minor,  but  quite  important  difference  among 
methods  centers  on  the  explicit  definition,  or  not,  of  a  target. 
Our  proposed  algorithm  requires  definition  of  seed  and  target 
regions  to  define  connectivity.  In  this  sense,  connectivity  is 
defined  similarly  by  the  approach  of  Hageman  et  al.  [12]  and 
MC  approaches  in  which  the  target  is  defined.  Such  condition¬ 
ing  on  the  target  eliminates  distal  fading  from  the  seed  to  the 
target,  as  can  be  seen  from  the  track  densities  being  the  same 
at  the  seed  and  target  regions  (Fig.  4).  The  other  approaches 
discussed  above  do  not  require  an  explicit  target  and  therefore 
result  in  maps  of  connectivity  between  a  seed  region  and  all 
other  voxels  which  fades  with  distance  from  the  seed.  For 
example,  in  graph-based  methods,  the  probability  of  a  given 
connection  is  derived  from  the  product  of  edge  weights  along 
a  path.  As  each  edge  weight  is  less  than  or  equal  to  1,  the 
probability  decreases  with  path  length,  and  the  connectivity 
maps  exhibit  distal  fading.  The  graph-based  methods  will, 
however,  always  define  a  pathway  between  the  seed  region 
and  any  given  voxel  even  if  the  connectivity  is  vanishingly 
small. 

A  practical  distinction  is  in  terms  of  speed.  Fluid  flow 
approaches  may  have  computation  times  of  hours  [12]  while 
our  proposed  method  takes  only  seconds.  The  proposed 
method  also  compares  favorably  with  graph-based  methods, 
especially  on  larger  computation  volumes  (Fig.  A3).  The  need 
to  run  the  computation  over  numerous  seed-target  pairs 
tempers  the  speed  advantage  but  can  be  addressed  by  calculat¬ 
ing  the  result  of  each  pair  in  a  parallel  fashion. 

C.  Limitations  and  Uncertainties 

The  simulation  shown  in  Figs.  3  and  4  demonstrates  that 
discretization  error — due  to  a  combination  of  discretization  of 
the  fODF  and  restriction  of  results  to  the  discrete  imaging 
grid — induces  systematic  errors  in  the  estimation  of  track 
densities,  but  none  so  severe  as  to  cause  failure  to  map  a 
connection.  In  addition  to  this  systematic  error,  the  restriction 
of  results  to  a  discrete  voxel  grid  can  be  a  liability  when 
mapping  pathways  that  are  small  compared  to  the  voxel 
dimension.  As  seen  in  Fig.  A2,  pathways  that  are  small 
compared  to  the  voxel  size  are  not  as  conspicuous  as  longer 
pathways.  Sampling  at  higher  resolution  can  ameliorate  this 
problem  to  an  extent.  However,  signal  to  noise  ratio  considera¬ 
tions  limit  acquiring  data  at  spacing  finer  than  approximately  1 
mm  cubed.  Interpolation  may  also  be  used  but  has  not  been 
examined  with  regard  to  the  proposed  method.  Conversely,  the 
performance  of  the  proposed  method  is  particularly  appealing 
for  long  pathways.  Distal  fading  of  track  densities  can  lead  to 


excruciatingly  long  computation  times  with  MC  methods  for 
long  pathways. 

Another  limitation  related  to  presenting  results  in  discrete 
voxels  is  that  the  extent  of  the  track  density  map  does  not 
appear  as  sharp  as  what  one  may  expect  from  streamline 
tractography.  As  the  MC  method  shows  this  to  a  qualitatively 
similar  extent,  we  believe  this  appearance  stems  primarily 
from  the  use  of  voxel-wise  track  density  maps  rather  than 
issues  such  as  discretization  of  the  fODF  in  the  calculation  of 
the  PDE  results. 

Unfortunately,  our  method  does  not  include  an  automated 
means  of  setting  thresholds  for  optimal  display  of  the  track 
density  map.  The  difficulties  are  illustrated  in  Fig.  A4.  Long 
and  short  association  fibers  are  shown  at  different  thresholds. 
Figures  with  minimum  threshold  set  to  zero  (Fig.  A4(a),  (f)) 
show  that  PDE  tractography  always  finds  a  nonzero,  but  often 
small,  value  for  track  density.  Finding  a  specific  threshold  that 
eliminates  the  implausibly  low  connectivity  values  depends 
largely  on  the  extent  of  the  white  matter  pathway,  as  do  the 
differences  in  how  the  tract  sharpens  and  diminishes  with 
threshold.  The  optimal  threshold  range  for  the  long  fiber  is 
lower  than  that  for  the  short  fiber.  The  value  of  track  density 
in  any  given  voxel  is  the  fraction  of  tracks,  leaving  the  seed 
and  hitting  the  target  but  not  the  boundary,  that  pass  through 
that  voxel.  As  can  be  seen  from  Figs.  A4(e)  and  (j),  the  track 
density  is  necessarily  high  near  the  seed  and  target.  However, 
in  intervening  voxels,  the  track  density  will  depend  on  a 
number  of  factors  specific  to  the  pathway,  such  as  the  length 
and  dispersion  in  space.  The  difference  in  track  density  at  the 
ends  and  middle  of  a  pathway  can  also  be  seen  in  the  “U” 
shape  of  the  track  density  plotted  against  position  in  the 
single-fiber  phantom  (Fig.  4(a)).  The  track  density  at  the 
middle  of  the  pathway  is  especially  low  in  Fig.  4(a)  because 
the  entire  imaging  volume  is  used  as  the  mask,  providing  a 
large  volume  for  tracks  to  disperse.  If,  at  the  other  extreme, 
the  mask  included  only  the  fiber,  the  track  density  in  the 
middle  of  the  pathway  would  be  the  same  as  at  the  seed  and 
target. 

Another  limitation  of  this  approach  is  that  there  is  no 
apparent  way  to  impose  geometric  constraints  on  the  tracks. 
The  formalism  considers  the  net  number  of  tracks  passing 
through  a  voxel,  but  the  geometry  of  the  tracks  counted  may 
be  unphysical.  Impossibly  sharp  bends  and  multiple  visits  to 
the  same  voxel  are,  for  example,  not  forbidden.  An  adjustment 
of  the  solution  algorithm  could,  in  principle,  impose  some 
local  constraints.  For  example,  the  continuity  equation  (2) 
could  be  solved  by  finding  values  in  the  immediate  neighbor¬ 
hood  of  the  seed  (2)  and  iterating  until  a  solution  is  found  at 
all  other  voxels  in  an  approach  similar  front  propagation 
methods  proposed  by  Parker  et  al  [13].  Constraints  could  be 
imposed  by  imposing  a  cost  on  the  curvature  of  the  front, 
thereby  preventing  implausibly  sharp  turns.  The  simplest 
implementation  could,  for  example,  start  at  the  seed  and  solve 
for  the  track  densities  in  voxels  neighboring  the  seed  using 
continuity  (2).  We  could  then  proceed  to  the  neighbor  with 
highest  density,  solve  for  the  densities  in  its  neighbors  and 
repeat.  As  we  progress  from  voxel  to  voxel,  we  can  enforce 
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physical  constraints  such  as  disallowing  high  bending  angles. 
Once  a  stopping  criterion  is  met,  such  as  arrival  at  the  target, 
we  would  have  an  approximation  of  the  track  density,  similar 
to  a  first  iteration  of  the  Gauss-Seidel  method  [20].  As  the 
approach  is  inherently  iterative,  there  may  be  an  undesirable 
cost  in  terms  of  computational  burden. 

An  uncertainty  that  requires  further  investigation  is  the 
cause  of  the  difference  between  the  PDE  and  MC  results  when 
tracking  the  transcallosal  (Fig.  11)  and  long  association  (Fig. 
Al)  pathways.  For  the  transcallosal  pathway,  while  the  MC 
results  show  a  track  density  through  the  pons  comparable  in 
magnitude  to  the  track  density  through  the  corpus  callosum, 
the  PDE  results  show  a  stronger  connection  through  the  corpus 
callosum.  This  difference  in  the  PDE  results  is  not  simply  due 
to  image  thresholding  or  presentation.  While  the  ratio  of  track 
density  in  the  pons  to  that  in  the  corpus  callosum  is  ~1  in  the 
MC  result,  it  is  only  -0.01  in  the  PDE  result.  The  cause  of  the 
difference  is  unclear,  but  we  did  investigate  several  possibili¬ 
ties.  We  made  modifications  to  the  MC  method  to  mimic 
features  of  the  PDE  algorithm.  To  approximate  the  discretiza¬ 
tion  of  the  fODF  in  the  PDE  method,  we  effectively  smeared 
the  fODF  by  using  a  smaller  maximum  order  for  spherical 
harmonics.  We  eliminated  the  physical  constraints  in  the  MC 
routine  that  disallow  backtracking  (bending  angles  greater 
than  90°),  which  is  allowed  in  the  PDE  method.  We  ran  the 


MC  method  with  these  changes  separately  and  together,  but 
saw  no  decrease  in  the  track  density  through  the  pons.  The 
same  mask  was  used  for  both  methods.  A  similar  difference 
was  seen  between  the  PDE  and  MC  results  in  the  superior 
longitudinal  fasciculus.  Given  that  the  PDE  method  was 
designed  to  mimic  the  behavior  of  the  MC  method,  we  are 
puzzled  by  the  differences  seen  in  Figs.  11  and  Al.  Further 
work  will  examine  a  wider  range  of  pathways  and  unravel  the 
reasons  driving  differences  seen  between  MC  and  PDE 
tractography. 

VI.  Conclusion 

We  have  introduced  a  generic,  probabilistic  formalism  for 
tractography  that  inherently  accounts  for  global  information 
while  having  a  fast  implementation.  The  advantages,  particu¬ 
larly  with  regard  to  long  pathways  and  pathways  along  the 
subcortical  surface  are  demonstrated  in  simulations  and  in 
vivo.  The  results  are  largely  equivalent  to  a  standard  Monte 
Carlo  approach,  but  faster  by  orders  of  magnitude.  The 
approach  promises  to  enable  whole-brain  tractography  with 
little  computational  burden. 

Appendix 

We  present  Figs.  A1-A4  referenced  in  the  text. 


Fig.  Al.  Maximum  intensity  projection  of  the  superior  longitudinal  fasciculus  generated  by  PDE  (a-c),  Monte  Carlo  (d-f),  and  shortest  path  (g-i)  methods. 
Intensity  ranges  were  set  to  [0.06,  0.25]  (a-c)  and  [0.01,  1]  (d-f). 
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Fig.  A2.  Maximum  intensity  projection  of  a  short  association  fiber  (U-fiber)  generated  by  PDE  (a,  b),  Monte  Carlo  (c,  d),  and  shortest  path  (e,  f)  methods. 
Intensity  ranges  were  set  to  [0.05,  0.5]  (a,  b)  and  [0.4,  1]  (c,  d).  Since  the  length  of  the  fiber  is  short  compared  to  the  voxel  resolution,  we  see  that  the  features  of 
the  track  are  not  as  distinct  as  those  of  longer  fibers. 


Mask  Volume  (voxels) 

Fig.  A3.  Log-log  plot  of  PDE  (x)  and  shortest-path  (o)  computation  time  vs.  number  of  voxels  in  the  mask.  PDE  computation  time  varies  linearly  with 
computational  volume,  while  shortest-path  time  varies  quadratically.  Linear  and  quadratic  dependence  are  indicated  by  the  dashed  and  solid  lines,  respectively. 
Red  symbols  correspond  to  times  for  computational  volumes  approximately  equal  to  that  of  the  in  vivo  data. 


Fig.  A4.  Comparison  of  long  and  short  association  fibers  at  different  thresholds,  (a)-(e)  Maximum  intensity  projections  of  the  superior  longitudinal  fasciculus  at 
minimum  intensity  thresholds  of  0,  0.05,  0.1,  0.15,  and  0.2,  respectively.  Maximum  intensity  is  set  to  0.25  for  all  five  panels.  (f)-(j)  Maximum  intensity 
projections  of  a  short  U-fiber  at  thresholds  of  0,  0.1,  0.2,  0.3,  and  0.4,  respectively.  Maximum  intensity  is  set  to  0.5  for  all  five  panels. 
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Target  audience:  Diffusion  MRI,  tractography,  intracranial  electrodes/EEG,  epilepsy 

Purpose:  Diffusion  weighted  imaging  (DWI)  can  non-invasively  infer  neural  tractography  and  connectivity,  which  can  reasonably  delineate  major 
known  pathways.  Regarding  accuracy,  there  is  continued  interest  is  how  well  DWI  measures  of  imaging  connectivity  (IC)  correlate  with 
electrophysiological  connectivity  (EC) — a  question  of  clinical  relevance.  One  standard  for  measuring  EC  is  cortico-cortical  evoked  potentials 
(CCEPs),  which  uses  intracranial  electrodes  to  stimulate  and  record  electrical  activity1.  Here,  we  compare  IC  and  EC  between  pairs  of  regions  in  the 
brain  varying  in  space.  Although  significant  correlation  between  these  two  measures  has  been  reported2,  our  results  indicate  that  the  relation  between 
IC  and  EC  is  modest  and  still  uncertain. 

Methods:  CCEP  data  was  obtained  in  epilepsy  patients  from  both  surgically  implanted  deep  stereotactic  electroencephalography  (SEEG)  electrodes, 
and  subdural  grid  electrodes.  EC  was  evaluated  by  directly  stimulating  an  electrode  contact  pair  and  recording  voltage  at  all  other  contacts.  A 
measure  of  EC  between  contacts  X  and  Y  was  taken  to  be  the  root  mean  square  (RMS)  of  the  recordings  at  Y,  from  10  to  20  ms  after  X  was 
stimulated. 

Prior  to  implantation,  the  patient  underwent  pre-procedural  imaging  at  3T  (2.5mm  voxels),  which  included  61  direction  diffusion  imaging 
(HARDI)3,  and  resting-state  connectivity.  All  electrode  positions  were  coregistered  to  anatomic  imaging.  We  have  developed  a  fast  probabilistic 
tracking  method4,  which  was  used  to  compute  pathways  from  each  stimulus  contact  to  every  other  contact,  to  parallel  the  CCEPs  stimulus.  Various  IC 
measures  include  either  mean  tract  density  along  pathways,  mean  tract  FA  or  TD,  or  mean  component  along  FOD.  Data  has  been  obtained  from  four 
patients  so  far. 


Results:  Fig.  1  shows  a  scatter  plot  (log-log)  of  the  relation  of  EC  and  IC  (DWI  connectivity)  for  over  600  contacts  in  4  patients;  each  data  point 
represents  a  pair  of  points  in  the  brain.  Although  there  is  a  significant  correlation,  there  are  many  points  showing  high  IC  with  low  EC,  and  visa 
versa.  Fig.  2  shows  a  similar  scatter  plot  comparing  EC  with  resting-state  connectivity,  which  shows  only  a  mild  correlation  of  EC  with  IC.  Fig.3 
shows  a  significant  distance  effect  whereby  both  IC  and  EC  diminish  with  distance,  an  effect  that  further  confounds  the  direct  correlation  of  EC  with 
IC. 

Discussion:  Today,  IC  is  extensively  used  to  study  the  brain,  but  there  is  little  evidence  supporting  various  measures  of  IC.  We  attempt  to  derive  a 
measure  of  IC  by  comparing  values  against  the  presumed  gold-standard  of  EC.  We  find  only  a  modest  correlation  with  DWI  and  a  minimal 
correlation  with  rsMRI.  Many  pairs  of  points  in  the  brain  have  high  IC  and  low  EC,  or  have  low  IC  with  high  EC.  We  have  explored  many  variations 
of  both  EC  and  IC  measures,  and  this  trend  is  always  observed.  Possible  causes  include  incompatibly  low  MRI  resolution,  distance  effects  of  both  IC 
and  EC,  poor  co-localization  of  MRI  voxels  with  electrical  contact  points,  or  intrinsic  brain  functionality  whereby  a  solid  stmctural  path  between  two 
points  does  not  necessarily  imply  a  strong  concordant  electrophysical  response.  Our  results  lead  us  to  believe  that  the  connection  between  EC  and  IC 
is  tentative.  This  discordance  could  be  inherent  to  brain  physiology,  or  it  could  indicate  that  our  current  methodology  is  not  fully  understood. 
Conclusion:  An  extensive  comparison  of  diffusion  MRI  and  electrophysiological  data  indicates  an  uncertain  correlation  between  these  two  measures 
of  brain  connectivity. 
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